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Abstract. The Szekeres family of inhomogeneous solutions, which are defined by six arbitrary 
metric functions, offers a wide range of possibilities for modelling cosmic structure. Here we 
present a model construction procedure for the quasispherical case using given data at initial 
and final times. Of the six arbitrary metric functions, the three which are common to both 
Szekeres and Lemaitre-Tolman models are determined by the model construction procedure 
of Krasinski & Hellaby. For the remaining three functions, which are unique to Szekeres 
models, we derive exact analytic expressions in terms of more physically intuitive quantities 
- density profiles and dipole orientation angles. Using MATLAB, we implement the model 
construction procedure and simulate the time evolution. 



Contents 



1 Introduction 1 

2 Lemaitre-Tolman models 3 

2.1 The LT metric 4 

2.2 Model construction: initial and final density profiles 5 

2.2.1 Coordinate choice 5 

2.2.2 Locating the worldlines of constant M 6 

3 Szekeres models 7 

3.1 The Szekeres metric 7 

3.2 Quasi-spherical case 7 

3.2.1 Riemann projection 7 

3.2.2 The function E and spatial foliations 8 

3.2.3 E'/E dipole 8 

3.2.4 Shell separation 10 

3.3 Singularities 11 

3.4 Regularity conditions 11 

4 Towards a model construction procedure 12 

4.1 Obtaining the arbitrary functions 12 

4.1.1 Density extrema 13 

4.1.2 Solving for E'/E 13 

4.1.3 Solving for S, P and Q 14 

4.2 The deviation function and origin behaviour of E'/E 15 

4.3 Expressing p ma x in terms of p m in and plt 17 

4.4 Evolution considerations 18 

4.4.1 LT density 18 

4.4.2 Reconstructing p m in from (E' / E) max during evolution 18 

4.4.3 Approximating R1/R2 18 

4.5 The algorithm 19 

5 Numerical simulations 20 

5.1 Implementation 20 

5.2 Model profiles 21 

5.2.1 Run #1 22 

5.2.2 Run #2 22 

5.2.3 Run #3 22 

5.3 Results and discussion 23 

5.3.1 Run #1 23 

5.3.2 Run #2 23 

5.3.3 Run #3 23 

6 Conclusions 24 



A Details of the LT model construction 

A.l Finding / and tf, 

A.2 Solving ip(x) = 

A. 3 Limiting values at M = 

A. 4 Reconstructing model evolution 



31 

31 
37 
37 

38 



1 Introduction 

Einstein's theory of General Relativity (GR) is widely regarded within the scientific commu- 
nity to be the best current description of gravitational phenomena, having been confirmed 
by a wide range of observational tests [1-6]. While this is not to say that it is necessarily 
a complete or final theory of gravity, as is evident by the numerous attempts to modify it 
[7-12], but rather, it does a good job of explaining, and in some cases predicting, obser- 
vations rather elegantly. It is thus of great importance to fully understand it's behaviour, 
and one key method is to investigate exact solutions. Assumptions of isotropy together 
with the Copernican Principle, and more particularly the cosmological principal, point to 
the Freedmann-Lemaitre-Robertson- Walker (FLRW) geometry as being a good approxima- 
tion to the behaviour of the Universe on the largest scales, with linear perturbation theory 
providing the first order corrections. And, in the mainstream of modern cosmology, almost 
all observations are analysed and interpreted within this framework. However, there are 
certainly scales on which exact non-linear GR is needed, possibly large scales. 

When attempting to specify a realistic model of the Universe, it is important to be 
able to relate quantities that characterise the model to observables in the real Universe. 
The traditional approach when studying dynamics is to specify the initial conditions of the 
system, i.e. initial data at some initial time, say t\, and then evolve it forward in time to 
some later instant, say t^- In the case of cosmology one may like to take decoupling as 
t±, and the present as t2, so as to be able to compare one's model with observations. The 
trouble is that our observations of the last scattering surface can only weakly constrain some 
of the initial conditions - others must be evolved backwards, from observations at the current 
time back to t±, which is often not possible. What ends up happening is one assumes an 
ansatz for the metric functions and attempts to 'tweak' them to fit observations. In addition, 
the evolution is very sensitive to velocity fluctuations, so it is extremely difficult to get a 
reasonable evolution from a data set at a single time. The ability to specify a model from 
a combination of both initial and final data is thus of great utility when trying to construct 
realistic models from observable data. 

Perhaps the most popular inhomogeneous exact solution used in cosmology today is 
the Lemaitre-Tolman (LT) model. First published in 1933 by Lemaitre [13-15], and then by 
Tolman [16] and later popularised by Bondi [17], it has since found numerous applications 
in the fields of astrophysics and cosmology. It is a spherically symmetric non-static solution 
to the Einstein Field Equations (EFEs) with a dust source, and can be thought of as an 
assembly of concentric spherical mass shells, each with their own evolution. Inspired by the 
work of Bonnor [18], Krasinski & Hellaby published a series of papers [19-22] (Hereafter 
referred to as Paper I, II, III and IV respectively) in which they developed a procedure 
for constructing LT models from a combination of initial and final data. In Paper I, the 
authors proved that any two spherically symmetric density profiles specified on any two 
constant time slices can be joined by a LT evolution, and exact implicit formulae for the 
arbitrary functions that determine the resulting LT evolution were given. Paper II extended 
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the procedure to determining the LT model which evolves from an initial velocity distribution 
to a final state that is determined either by a density distribution or a velocity distribution. 
Numerical examples of the evolution of various structures were investigated, and it was found 
that initial velocity fluctuations have a much bigger effect on the subsequent evolution than 
initial density fluctuations. Paper III applied the methods developed in previous papers to 
the formation of a galaxy with a central black hole. For this investigation, the first of its kind, 
the tools of GR were essential since it is not possible to describe a black hole in Newtonian 
gravity or as a perturbation around a FLRW background. Paper IV went on to generalise 
the model construction methods given in the previous three papers, to find the LT model 
which satisfies any two of the following 'boundary conditions': a simultaneous big bang, a 
homogeneous density or velocity distribution in the asymptotic future, a simultaneous big 
crunch, a simultaneous time of maximal expansion, a chosen density or velocity distribution in 
the asymptotic future, only growing or only decaying modes. Some of these new specification 
methods were used in [23] to model the Shapely Concentration and the Great Attractor. The 
importance of velocity profiles was again highlighted in [24] where it was demonstrated that 
an initial over-density can evolve into a void, given a suitable initial velocity profile. 

Another interesting family of inhomogeneous exact solutions are those found by Szekeres 
[25], in 1975. In general, these models have no symmetries (i.e. no killing- vectors [26]) and 
are defined by six arbitrary metric functions - representing a freedom to rescale the 'radial' 
coordinate and five degrees of freedom to model inhomogeneity. They are perhaps the most 
sophisticated exact solutions with a dust source, and offer exciting prospects for modelling 
fairly complex cosmic structures. There are two classes of Szekeres models, the LT-type 
(/3, z 7^ or Class I) and the Kantowski-Sachs (KS) type ((3, z = or Class II) 1 Bonnor [29, 30] 
showed an interior region of LT-type quasispherical Szekeres spacetime can be matched to 
the exterior Schwarzschild solution, even though the interior metric has no symmetry. Since 
the Schwarzschild solution does not contain any gravitational radiation, this implies that 
such Szekeres models do not radiate, and consequently proves the existence of configurations 
of collapsing dust clouds that have no symmetry and do not produce gravitational waves. 
Goode & Wainwright [31, 32] introduced a different representation of the Szekeres solutions 
in which many properties of both subfamilies can be considered together 2 . Furthermore, 
this formulation facilitates the separation of 'exact perturbations' from background FLRW 
dynamics. Recently, it has been used by Ishak & Peel [34] to study the evolution of large 
scale structure. Meures &: Bruni [35] recently considered the KS-Type Szekeres solutions 
with A / 0, originally obtained by Barrow &: Stein-Schabes [36], to model an arbitrary initial 
matter distribution along one line of sight. They re-parametrised the solutions into the Goode 
&; Wainwright representation, and gave exact solutions for the growing and decaying modes 
of the metric perturbation, assuming a flat ACDM background. 

Within the LT-type class there are three further subclasses, the most popular of which is 
the 'quasispherical' model, which we focus on here. The geometry of the 'quasi-pseudospherical' 
and 'quasiplanar' models is still poorly understood [37], perhaps due to our lack of under- 
standing of non-spherical gravity, and so, have not been explored for cosmological applica- 
tions. The quasispherical model can be thought of as a sequence of non-concentric mass 
shells, each with a density dipole distribution and its own evolution - a generalization of the 
LT model [38]. Soon after his publication of these solutions, Szekeres used the quasispherical 

x For a comprehensive review of the two classes of Szekeres models, see §19.6 in [27] or, for a more historical 
account, see §2.4 in [28]. 

2 The KS-Type was later shown to be a regular limit of the LT-Type [33] 
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model to study non-spherical gravitational collapse [39], and since then they have found cos- 
mological application in the study of light propagation [40] as well as structure formation [41] . 
Bolejko [42] used the quasispherical model to study the evolution and interaction of a void 
with an adjoining galaxy cluster. He found that small voids surrounded by large overdensi- 
ties evolve much slower than large isolated voids do. And similarly, large voids enhance the 
evolution of adjacent galaxy superclusters, causing them to evolve much faster than isolated 
ones. Sussman & Bolejko [43] presented an approach to describing the dynamics of Szekeres 
models in terms of 'quasi-local' scalar variables. In this formulation the field equations and 
basic physical and geometric quantities are formally identical to their corresponding expres- 
sions in the LT model, thus potentially allowing the generalisation of rigorous LT results 
to the non-spherical Szekeres geometry. They then used this formalism to investigate small 
dipole perturbations away from spherical symmetry, showing that such a configuration is in 
fact stable. In [44] Bolejko &; Sussman used the quasispherical model to construct several 
elongated supercluster-like structures with underdense regions between them, thus providing 
a reasonable coarse-grained description of present day cosmic structures. After averaging, 
the authors found that such a density distribution produced a spherical void profile which is 
roughly consistent with observations, without the need for dark energy. Also, by considering 
a non-spherical inhomogeneity, the definition of a "centre" becomes much more nuanced, and 
thus constraints placed by fitting observations on our position relative to this location become 
less restrictive. Mishra et al [45] calculate the redshift drift equation for an axially symmetric 
quasispherical Szekeres Swiss-cheese model, along the axis of symmetry. They compare these 
results to the drift in the ACDM model and some LT models, finding that they are a good 
discriminator between them. They go on to propose a method to fully constrain all degrees 
of freedom in an axially symmetric quasispherical Szekeres model from observable quantities. 
Apparent horizons for the quasispherical Szekeres model have recently been calculated by 
Krasinski & Bolejko [46]. 

The aim of this paper is to develop an algorithm by which one can construct realistic 
Szekeres models from some given data on the initial and final time surfaces. This will entail 
deriving expressions for the six arbitrary metric functions, which will completely define the 
Szekeres model between the initial and final time, from some physically more intuitive quan- 
tities. Here we choose to work with initial and final density data, although it is foreseeable 
that this approach can be extended to include velocity data. The structure of this paper is as 
follows: §2 presents the LT model and a summary of the construction procedure; §3 presents 
the Szekeres model; §4 presents the derivation of expressions for the arbitrary functions, and 
outlines the model construction procedure; §5 presents the results of the numerical simula- 
tions using the procedure outlined in the previous chapter. Conclusions are then presented 
in §6, along with some comments on future work. 

2 Lemaitre-Tolman models 

A brief review of the LT model is necessary because (a) many of its functions carry over to 
the Szekeres models unchanged, and (b) we are generalising a model construction method 
given for the LT model. 
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2.1 The LT metric 

This metric, in geometric units {G = c = 1), is 

T>/2 

ds 2 = -dt 2 + Y^~f dr2 + R2dn2 t 2 ' 1 ) 

where d£l 2 = d6 2 + sin 2 6d(f> 2 is the metric of a unit 2-sphere, and ' = J^. The function 
R = R(t,r) is known as the areal radius as it is related to the area of constant- (i, r) 2- 
surfaces, while / = f(r) is an arbitrary function that determines the type of evolution and 
the local geometry. The dust source is described by an energy-momentum tensor for a 
pressure-free perfect fluid, 

T ab = pu a u\ (2.2) 
that is comoving with the coordinate system, such that 

u a = 5?. (2.3) 
Applying the EFEs to the metric yields two expressions - an equation of motion, 

R * =f+ ™ + l R i, (2 . 4 ) 
and an expression for the energy density, 

2M' . . 

Kp= mi" (2 - 5) 

where ' = and M = M(r) is another arbitrary function that gives the gravitational mass 
within a comoving shell of 'radius' r. One can see that (2.4) is simply the Friedmann equa- 
tion for dust, except that / and M are functions of r. It is evident from (2.4) that f(r) also 
represents twice the local energy density per unit mass of the dust particles, so it is often 
written f(r) = 2E(r). 

The evolution of R depends on the value of /. When A = 0, the solutions of (2.4), in terms 
of a parameter rj, are 3 



Hyperbolic (/ > 0) 



M , 
R = — (cosh^ — 1) 

(smhr,- V ) = f3/2 ^ tb) (2.6) 



Parabolic (/ = 0) 



R = M— 
2 

V 3 t-t b 



(2.7) 



6 M 

3 Near the origin, where / — > and M — > 0, the evolution type is determined by the sign of Rf /M or 
f/M 2/3 



-4- 



Elliptic (/ < 0) 



M 

r = rr^yi 1 - cosr ?) 

(r? _ sin?? ) = ^l_A H (2.8) 

where tfe = tb(r) is the last arbitrary function - the 'bang time', which gives the local time 
of the initial singularity t = tj, (i.e. when R = on each worldline). So particle worldlines 
can emerge from the bang at different times, typically outer spheres before inner ones. Since 
/ = /( r )> it is entirely possible to have adjacent regions of hyperbolic and elliptic evolution. 
These regions will be connected by a parabolic shell (or extended region) at the boundary, 
since / is required to be continuous. A nice example of adjacent elliptic and hyperbolic re- 
gions is a re-collapsing dust cloud surrounded by an ever-expanding universe [47] . The time 
reversed parabolic and hyperbolic cases, obtained by writing (tb — t) instead of (t — tb), are 
also valid solutions. 



A useful expression for R' is [48] 

This is valid for all three expressions (2.6), (2.7) 4 , and (2.8). 

2.2 Model construction: initial and final density profiles 

The most obvious way to construct a LT model is to specify the arbitrary metric functions, 
M(r), f(r) and tb(r). Alternatively, one makes a coordinate choice and specifies the initial 
conditions - for example; the density p(t\,M) and velocity R(t\,M) distributions at some 
initial time t = t\. Since it is not always obvious what density distribution and model 
evolution will result from a particular choice of initial conditions, there are situations were it 
is preferable to determine the metric functions from a combination of initial and final data. In 
[19] the authors demonstrated that any two spherically symmetric density profiles specified 
on any two constant time slices can be joined by a LT evolution, and gave exact implicit 
formulas for the arbitrary functions that define the resulting LT model. Here we present a 
slight modification to this procedure, which will be useful for the proposed Szekeres model 
construction procedure presented in §4. 

2.2.1 Coordinate choice 

The original formulation of [19] made the coordinate choice r = M(r), however, here we find 
it convenient to make the choice 

r = R(t 2 ,r), (2.10) 

and thus 

R'\t=t 2 = I- (2.11) 

4 In (2.9), one need not set /'// = for the parabolic case, as claimed in [49]. On the boundary between 
hyperbolic and elliptical regions, one has / = and /' 7^ 0. An /' term remains in (2.9) if the parabolic limit 
is taken correctly. 



t' h + 



M 3f 



t b ) 



R. 



(2.9) 
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There are a couple of reasons why this choice is preferable. Firstly, the origin limit calculations 
for E' /E, shown in §4.2, are simpler when R' = 1. Secondly, the choice is less restrictive 
on the allowable density profiles which produce a finite origin value of E' /E 5 . And thirdly, 
r = M does not allow vacuum regions, M = constant, where the r coordinate would be 
degenerate. Choice (2.10) allows one to write the terms with 'radial' derivatives instead with 
respect to R 2 , as 

d d „. dR %r , dM . . 

5? - m * R «v = m (2 - 12) 

2.2.2 Locating the worldlines of constant M 

Our procedure for determining the LT metric functions differs from [19] only in how one finds 
corresponding values of R and M from the initial and final density profiles (essentially, how 
one determines M(r) with the choice (2.10)). Suppose the density distributions at the initial 
instant, t = t\, and the final instant, t = t2, are given by 

piGRO = p(h, Rt), p 2 (R 2 ) = P (t 2 , R 2 ), (2.13) 

where 



Ri{r) := R(t h r), i = l,2. (2.14) 

Since M is constant along particle worldlines, corresponding values of R (and p) at t\ and t 2 
can be found by rearranging (2.5) and integrating over R. One finds 

M(U, Ri) = M min + £ / 1 Pi R* 2 dR* i = 1, 2, (2.15) 

" J H min 

though M m i n and Rimin will typically be zero, unless there is a central black hole. For 
definiteness, it is assumed in the following that the final instant is later than the initial 
instant, and that the final density is smaller than the initial density, at the same M. That is 6 

p{t 2 ,M) < p{t u M), t 2 >h (2.16) 

This implies that matter has expanded along every world line, although the analysis can 
be easily adapted to the collapse situation. As a result of the assumption (2.16) one finds 
R 2 {M) > Ri(M) (see Eqn (3.2) of [19]). Thus, evaluating (2.15) with (2.13) allows one 
to determine the corresponding values of Ri and R 2 (see §4.4.3 for details on the ranges of 
integration). Once these corresponding values are determined, the procedure for finding / 
and tb follows exactly as in [19]. Since we do not wish to duplicate already published work, 
we refer the reader to the original article for details of the derivation. In summary, there are 
three main cases and two borderlines which can evolve between t\ and t 2 . The remainder 
of the procedure, as adapted to the above coordinate choice, is summarised in appendix A. 
This procedure is central to the Szekeres version that follows. 

5 The reason for this will become clear in §4.1.2, where we derive an expression for E' /E in terms of the 
density profiles 

6 Here p(t 2 , M) is used for p(t 2 , R 2 (r{M))) etc. 
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3 Szekeres models 



3.1 The Szekeres metric 

The LT-type Szekeres metric [25, 39] is 

e + / E L 

where e = ±1,0. For e + 1, the functions f(r) and R(t,r) have the same significance as in 
the LT model. The function E(r,p,q) may be written 



. S 
E(r,p,q) = - 



S \ S 



(3.2) 



where 5 = S(r), P = P(r) and Q = Q(r) are arbitrary functions, and S, P and Q have 
natural interpretations in the Riemann projection (see §3.2.1). The dust source is described 
by the energy-momentum tensor for a pressure-free perfect fluid that is comoving with the 
coordinates, such that 

T ab = pu a u b , u a = 5f. (3.3) 

Applying the EFEs to the metric yields two expressions. An equation of motion identical to 
(2.4), which has solutions (2.6) (2.7) (2.8), and an expression for the energy density 



Kp 



M' - 3M [ % 



m [sf - r (§)] 



(3.4) 



With e = +1, the function M(r) has the same interpretation as in the LT model - it gives 
the gravitational mass within a comoving shell of 'radius' r. Any equations from §2 involving 
only R, M, f, a and t, and their derivatives, continue to hold. 

3.2 Quasi-spherical case 

From the family of Szekeres models it is the e = +1 quasi-spherical case that has received the 
most attention in the field of cosmology. They have found applications in the study of the 
early Universe [50, 51], structure formation [41, 42], the dimming of the supernovae [52, 53], 
light propagation [40], CMB observations [54] and volume averaging [55]. The e < cases 
have been far less investigated [37] , the geometry is still not well understood, and they are yet 
to find cosmological application. While these models may prove of some use in future, in this 
paper we focus attention on the quasi-spherical case, in which we expect the generalisation 
of the LT model construction procedure to be easiest. 

3.2.1 Riemann projection 

In (3.1) the metric component (dp 2 + dq 2 )/E 2 is a unit 2-sphere, plane or pseudo-2-sphere 
in Riemann projection. When e = +1, the p-q 2-surfaces are related to 6-(p 2-surfaces by one 
of the following transformations. Either 

cot (!) cos ((f)) ( ^— —) =cot ) sva(4>) (3.5) 



S 12/ ^ \ S 12 
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tan ( — J cos(0) \ \ = tan ( - ) sin(0), (3-6) 



or 



s J V 2 / V s J \2 / 

with similar expressions [37] for e < 0. The transformed 2-metric is then 

ds 2 = R 2 (d9 2 + sm 2 {6)dct) 2 ). (3.7) 

Each of the spherical transformations, (3.5) and (3.6), cover the entire p-q plane with < 9 < 
7T and < (p < 2tt, so either can be used. See [37] for illustrations of the various projections. 

3.2.2 The function E and spatial foliations 

Since the transformations (3.5) - (3.6) are ill-defined if S = 0, it is normal to assume S > 0. 
Then it is clear from (3.2) that when e = +1 we have E > 0. Differentiating (3.2) with 
respect to r and setting E' = gives 



Q \ 1 2 

P- It- 
s' 



+ 



Q-Q'|) 



l2 = W^W +£l . (3. 8 ) 



This is the equation of a circle in the p-q plane centred at the point (p, q) = (P — P'S/S', Q — 
Q'S/S 1 ) with radius 5 V / (P' 2 + Q l2 )/S' 2 + e. When S' > 0, one finds E' < inside the 
circle and E' > outside [38]. The locus (3.8) has the effect of creating poles and zeros 
in the function E'/E, which appears in g rr and p, see (3.1) and (3.4), and thus plays an 
important role in the geometry and physics of the model. The constant-(t, r) 2-surfaces are 
multiplied by the factor R = R(t,r), and hence the r-p-q 3-spaces are constructed from a 
sequence of constant-(i, r) 2-spheres which have a radius R. The g rr component of the metric 
(3.1) is affected by p-q variations via the function E'/E (3.11), which means that the radial 
separation between two neighbouring spheres of constant-(t, r) also has p-q variation, and 
thus these 2-spheres are interpreted as being arranged non-symmetrically relative to each 
other. 

3.2.3 E'/E dipole 

Here we show that the function E'/E has a dipole variation around each constant- (t, r) 2- 
sphere, with the extrema located at antipodal points, and E'/E = on the 'equator' between 
the two 'poles'. In order to illustrate the variation of E'/E over a constant-(i, r) 2-sphere 
(9, (^-coordinates are convenient. Applying either of the spherical transformations (3.5) or 
(3.6) to (3.2) and its derivative, yields 

E=-—^-:, (3.9) 
1 — cos 9 



and 



The locus E'/E = is 



. S cos 9 + sin#(P cos<^ + Q sin 

E = 

1 — cos 9 

E' S' cos 9 + sin 9 (P' cos <f> + Q' sin 



1-cosO (3 ' 10) 



E S < 311) 



S' cos 9 + P' sin 9 cos dp + Q' sin 9 sin (j) = (3. 12) 
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and applying the rectangular transformations 



x = sin#sin0, y = sin#cos</> z = cos9 (3.13) 
allows for a natural interpretation 

P'x + Q'y + S'z = (3.14) 

which is the equation of an arbitrary plane passing through (0, 0, 0). The locus E'/E = is 
then the intersection of this plane with the unit 2-sphere, which is a great circle. The unit 
normal to the plane (3.14) is 

(315) 

By setting (3.11) equal to a constant and employing the rectangular transformations men- 
tioned above, one finds the loci E' /E =constant 

P'x + Q'y + S'z = kS k = constant (3.16) 

to be the equation of arbitrary planes parallel to (3.14). This implies that all the loci 
E'/E =constant are small circles parallel to the E' = great circle (3.12). The location 
on the 2-sphere of the E'/E extrema (i.e. the poles) are found by finding where the partial 
derivatives of E'/E, with respect to 6 and <fi, are equal to zero. Denoting the location of the 
extrema by (9 e ,(j> e ) one finds 

d{E'/E) _ sm9(P' sine/) - Q' cos <j>) 
bM> ~ S ~ 

^tan</> e = | 7 (3.17) 
P' 

=^cos</> e = ei — =, ei = ±1 (3.18) 

d(E'/E) _ S' sin 8 - P' cos Ocoscp- Q' cos 9 sin 
89 ~ ^ S 

cos 6 e = e 2 ; e 2 = ±1 (3.20) 

\/S' 2 + P' 2 + Q' 2 

Applying the transformation (3.13) to the expressions found above gives the location of the 
extrema in rectangular coordinates. One finds 

(x e ,y e , z e ) = €2 — , (3.21) 

v ' W( p ') 2 + (Q') 2 + (S') 2 

This vector points in the same direction as the vector normal to the plane of the great circle 
(3.15), so the extrema of E'/E are located at the poles of E'/E = great circle. The extreme 
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values of E'/E are found by substituting the location of the extrema (3.18) (3.20) into (3.11), 
which gives 

(e'\ Vs' 2 + P' 2 + Q' 2 



I -p 1 = s — (3-22) 

\ / extreme 

The parameters €2 and ei are not independent 7 . The relationship between the two can be 
found by relating the sign of (3.19) to that of (3.20), as follows. Noting that 9 is defined on 
the interval [0, tt], we have sine? > 0, and hence 

sign(cos# e ) = sign(tan6* e ) 
=> sign(e 2 S') = sign(ei/S') 

ei = e 2 . (3.23) 

This seems reasonable since we expect only two extrema. From (3.22) and (3.23), and 
assuming that 5 > always holds, we deduce that e% = e 2 = 1 corresponds to a dipole 
minimum, while the dipole maximum is given by e% = e 2 = — 1- From this, it is clear that 

/ max V / min 

and at the dipole maximum, the E'/E value is then 



while the orientation angles are 



P' 

COS (bmax = = (3.26) 

a/P' 2 + Q' 2 
S" 

cos^ max = , (3.27) 

Sketches of a dipole and the variation of its orientation can be found in section 4.4 of [47]. 
3.2.4 Shell separation 

As pointed out in §3.2.2, the g rr component of the metric (3.1) is sensitive to p-q (and hence 
9-(j>) variations via the function E'/E, and hence the const ant- (f, r) 2-surfaces are interpreted 
as being arranged non-symmetrically. Writing the radial separation as 

R/ _ R E[ = Rf + R S' cos 9 + P' sin cos <f> + Q' sin 6 sin 

we see that RE' / E is a correction term to the spherically symmetric radial separation, 
R', that an LT model would have. The minimum radial separation between neighbouring 
constant-r shells obviously occurs where E'/E is maximum. Hence, e = +1 Szekeres 3-spaces 
are interpreted as being constructed from a sequence of non-concentric 2-spheres, with each 
shell having the exact density distribution required to generate a spherical field around the 
new centre. The dipole variation in E'/E around each constant-(t, r) 2-sphere causes a dipole 



7 The presentation in [38] misses this point 
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variation in the density (3.4). At the equator, where E' /E = 0, we see that (3.4) reduces to 
the form 

2M' 

k PLT = ^ (3.29) 

which is identical to the expression for density in LT models (2.5). In the context of Szekeres 
models, we refer to the equatorial density (3.29) as the LT-density, denoted by the subscript 
LT. 

3.3 Singularities 

Szekeres models contain the same singularities as in LT - those of the bang, the crunch and 
shell crossings. The Kretschmann scalar is 

K = k 2 (i& v - ^ PAV p + 3p 2 ^j + (2A + np) (3.30) 

where 

6M 

WAV = -pg" (3.31) 

is like a "mean density" 8 interior to the comoving shell of constant-r. The dynamics of R 
in the Szekeres model is identical to that of the LT model (2.4), and hence the bang and 
crunch singularities are the same. Shell crossings occur when inner shells of matter pass outer 
ones, causing the density to diverge and the radial coordinate, r, to become degenerate. In 
Szekeres models they are more complicated than in the spherically symmetric LT case, as non- 
concentricities cause constant-r shells to pass through each other gradually, and one shell may 
intersect many others at any given time. So to avoid shell crossing one requires, in addition to 
the LT conditions, further constraints on the radial derivatives of the arbitrary functions S' , 
P 1 and Q' . Necessary and sufficient conditions to completely avoid shell crossings in Szekeres 
models are given in [38] for the A = case. 

3.4 Regularity conditions 

For the metric (3.1) to retain Lorentzian signature ( h ++), the g rr component must 

always remain positive, and thus e + / > is required, with the equality only occurring 
where R' — RE' / E = 0. In the e = +1 case the origin is the locus r = ro where the 2-sphere 
radius vanishes, i.e. R(t, ro) = V t so that R(t, ro) = 0, R(t, ro) = etc. This is not a 
centre of symmetry, but a locus where the shear and electric Weyl tensors vanish (assuming 
a regular origin). Regular origins require that on any constant t surface away from the 
bang or crunch, the density (3.4) and curvature (3.30) remain finite, and the time evolution 
at r = ro should be a smooth continuation of the immediate neighbourhood. Conditions 
on the arbitrary functions were investigated in [38] by evaluating the limit r — > ro of the 
aforementioned quantities and ensuring they are well behaved. They found the regularity 
conditions require that near the origin 

M ~ R 3 f ~ R 2 

S ~ R n P ~ R n Q~ R n < n < 1 (3.32) 

Of the five conditions above, the first two are the same as in LT models. Spatial extrema in R 
may occur in the constant t 3-spaces of models with certain topologies. For example, closed 



see e.g. [33], also c.f. the quasi-local average of [43] 
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spatial sections can have R increasing away from the origin until some point, say r = r m , 
where R is maximal and beyond which R is decreasing toward a second origin. The conditions 
to ensure regularity of these loci were investigated in [38] . Ensuring no shell crossings at the 
maximum, R' = requires 

M'{t, r m ) = f'(t, r m ) = t' B (t, r m ) = 

S'(t, r m ) = P'(t, r m ) = Q'(t, r m ) = V t (3.33) 

Furthermore, to avoid any surface layers at r = r m requires 

f = -e (3.34) 

With these conditions met, the density (3.4) and the g rr component of the metric (3.1) remain 
positive and finite, thus ensuring regular extrema. 

4 Towards a model construction procedure 

In order to construct a realistic Szekeres model one must specify the arbitrary functions 
associated with the model from some physical quantities. These functions, of which there are 
six in total, allow for a rescaling of the r coordinate, f = f(r), plus five physical degrees of 
freedom to model inhomogeneity. Once one has obtained expressions for these six functions 
the Szekeres model will be completely specified. Special attention must be paid to the origin 
as certain variables have a value of zero there. 

4.1 Obtaining the arbitrary functions 

Since the expression for the Szekeres equatorial density (3.29) corresponds exactly with the 
LT density expression (2.5), and the dynamics of R are identical in both Szekeres and LT 
models (2.4), the LT model construction procedure of Krasinski & Hellaby can be used 
to determine the arbitrary functions which are common in both models. Thus, specifying 
a Szekeres equatorial density profile (hereafter the "LT-density" ) at some initial and final 
time, t\ and t 2 , is sufficient to determine the arbitrary functions M, f and tf, according to the 
procedure described in §2.2. Note however that the radial coordinate choice is not as in [19]. 
The arbitrary functions S, P and Q require further knowledge of the intensity and orientation 
of the dipole. These dipole parameters can be specified on only one of the 3-surfaces, t± or 
t 2 - we choose to do it at the later time as we expect one would know more about density 
detail at later times. Defining the density extrema on a particular constant r shell at time 
t 2 to be 

Pm.in{r) = min [p(t 2 ,r, 9, 4>)] (4.1) 
6,tf> 

and 

Pmaxij) = max [p(t 2 ,r, 0,<p)] (4.2) 

0,0 

we specify the following: 

• Plt,i(Ri) LT-density profile at t = ti 

• pLTpiRi) LT-density profile at t = t 2 
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• Pmin{R2) Density minimum profile at t = t-2 

• Pmin {R2) Density minimum orientation angle, 9, at t = t 2 

• ^pmini-R-i) Density minimum orientation angle, <fi, at t = ti 

From these quantities we will extract expressions for the arbitrary metric functions, S, P 
and Q. 

4.1.1 Density extrema 

Previous literature ([38] - Equation 69) claimed to show the derivative of the density with 
respect to E'/E to be negative (i.e. p x < with x = E'/E), implying that E'/E is a 
minimum at a density maximum, and a maximum at a density minimum. Moreover, since 
E'/E is a dipole (see §3.2.3), it implies that a density maximum corresponds to a negative 
E'/E value, and a density minimum to a positive E'/E (i.e. E'/E\ Pmax = E'/E\ m i n < 
and E' /E\ Pmin = E' /E\ max > 0). This is in fact not the case. Taking the derivative of the 
density (3.4) with respect to x = E'/E, and using (3.29) and (3.31), we find 

{plt ~ PAv) 



RR' 



(4.3) 



{R> - Rx) 2 

Clearly the sign of p, x is not always negative, as was previously thought, but rather it depends 
on the sign of R'{plt — PAv), since R/{R' — Rx) 2 > always holds. Thus, E'/E\ max only 
occurs at the density minimum if R'(plt — PAv) < 0. 

4.1.2 Solving for E'/E 

In order to express E'/E in terms of density parameters, we begin by rearranging (3.4) to 
give 

E'\ k P R 2 R'-2M' 



(4.4) 

E J upR* -6M y ; 

Then, dividing both numerator and denominator by kR? and substituting in (3.31) and 
(3.29), leads to 

'E'\ R' (p-plt 



E J R \p-p AV J ' (45) 

The anti-symmetric property of E'/E (3.24) allows one to express E'/E\ max in terms of 
E'/E\ Pmin , regardless of the sign of R'(plt — PAv)- Substituting (4.5), with p = p m i n , into 
(3.24) gives 

—\ = — ( Pmin - PLT ^ (4 6) 

E ) max R \Pmin ~ PAV J 

However, when relating the orientation angles of the density-dipole to the E'/E dipole, it is 
necessary to track the sign of R'(plt — PAv)- Since the poles of E'/E (and p) are antipodes 
on the 2-sphere, the orientation angles are easily related. 

R'ipLT - PAV)>0: Qmax = V - B Pmm , (j) max =7T + <j) pmin (4.7) 

R'ipLT - PAV) < : Qmax = Pmin , <j) max =4> Pmin - (4.8) 

We now have expressions for E'/E\ max , max , <pmax as functions of v in terms of the physical 
quantities plt, PAV, Pmin, Pmin and (j) Pmin at time t 2 . We can thus go on to define the 
arbitrary functions, S, P and Q, using these expressions. 
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4.1.3 Solving for S, P and Q 

Equations (3.26), (3.27) and (3.25) suggest that the radial derivatives of the three arbitrary 
functions, S', P' and Q', can be solved for in terms of the dipole orientation angles, 9 max 
and (pmax, and E'/E max . The metric functions, S P and Q, can then found by integrating 
the expressions for S' , P' and Q' over a suitable radial coordinate. Solving the system of 
equations (3.25) (3.26) (3.27), we find the following. Multiplying (3.25) by S, squaring both 
sides and rearranging leads to 



incur 



and substituting this expression (4.9) into (3.27), and using (3.25), eliminates the S' terms, 
to give 



•&y s 2 -p' 2 -q 

, hi I max ^ 



COS0 max = — I -gTT . (4.10) 

V E J max 

Multiplying (4.10) by the RHS denominator, squaring and rearranging yields an expression 
for Q' 2 in terms of only P' and S, 

Q' 2 =(^) S 2 sin9 2 nax -P' 2 . (4.11) 

V / max 

Now Q' 2 can be entirely eliminated from (3.26) by substituting it into (4.11), and rearranging 
the result then gives an expression for P', 

P' = - sin 6 max cos <j) max (^) S. (4.12) 

\ / max 

By substituting (4.12) into (3.26) one eliminates P', and solving for Q' one finds 

Q' = 6 3 sin0 max .sin ( /w(^) S e 3 = ±1 (4.13) 

V / max 

From (3.25) and (3.27) it is easily seen that 

S' = -cos9 max (?p) S (4.14) 

\ / max 

In (3.25) - (3.27) the Q' terms are all squared and as a result, the sign of Q' in (4.13) is 
undetermined. To lift this degeneracy one uses (3.17) (4.12) and (4.13), obtaining 

f , _Q' 

tan (pmax — ~pj 

= -e 3 tan<j) max , (4.15) 

and it is clear from (4.15) that e 3 = —1. Hence the final equations for S' , P' and Q' are given 
by 

S > = - C o S ma J S (4.16) 

V / max 

P' = -coscj) max sin9 max (—J S (4.17) 

V / max 

Q' = - sin (j) max sin 9 max (^) 5 (4.18) 

\ / max 



- 14 - 



Expressions for the three arbitrary functions, S(r) P{r) and Q(r), are now easily obtained 
by integrating (4.16), (4.17) and (4.18) over the radial coordinate. Since all the quantities to 
be integrated are specified in terms of R 2 , it is sensible to integrate these expressions over 
the same radial coordinate. We find 



S(R 2 ) = So exp 



«2 / E > 

cos max ( — ) dR% 



E , 

' max 



S(R mm ) = S (4.19) 



fR2 / T?l s 

P(R 2 ) = p Q - j cos (j) max sin 6 max I — ) SdR* 2 (4.20) 



rn / n 



R2 ^ ( E > 



Q{R2) = Qq-J sin (p max sin 6 max ( — j SdR* 2 (4.21) 

J Rmin \ ' max 

where Pq and Q$ are initial values which can be set to zero. Since a constant rescaling of S 
has no physical effect on the metric, one is free to choose the value of Sq. For simplicity, we 
will use Sq = 1. Having obtained analytic expressions for the metric functions in terms of 
physical quantities, one can completely determine the Szekeres metric. While the quantities 
Smax, 4>max and {E' / E) max are not physical, they are related to the physical quantities 9 Pmin , 
4> P mi n i Pmin, PLT and p AV , as was shown in §4.1.2. 

4.2 The deviation function and origin behaviour of E' / E 

As the origin is approached, by l'Hopital's rule pav Plt, and we also expect p m in — > Plt 
since a pointlike dipole seems unphysical 9 . Thus, it is apparent from (4.6) that the origin 
value of (E'/E) max goes to 0/0. The conditions that ensure this quantity is finite at the 
origin are found by writing the density profiles, pur-, Pmin and pav, as series expansions 
about r = 0. On a constant time slice we make the co-ordinate choice r = R, and hence 
R' = 1. The series expansions for plt and p m in are then given by 

PLT = P0+PlR + P2R 2 + P3R 3 + - (4.22) 

p min = Co + CiR + C2R 2 + ( 3 R 3 + ... (4.23) 
Rearranging (3.29) and inserting (4.22) gives an expression for M', 

M' = ^[p R 2 + Pl R 3 + P2 R A + p 3 R 5 + ...], (4.24) 

which can be integrated to yield an expression for M, 

M = ^[p R 3 + -p^ + \p 2 tf> + IpsB? + ...], (4.25) 
6 4 5 6 

and substituting (4.25) into (3.31) gives 

3 3 3 

PAV = Po + -.PlR + TP2R 2 + t;Pz r3 + ■- ( 4 - 26 ) 

4 5 6 



9 Mathematically we merely require the terms R' - RE' /E = R'(l - RE'/R'E) and M' - 3ME'/E = 
M'(l - 3ME'/M'E) in (3.1) and (3.4) to be well behaved at the origin, so RE'/R'E -> is sufficient, since 
we already have M ~ R 3 there. 
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Now, inserting the series expansions (4.22), (4.23) and (4.26) into the expression for {E 1 /E) r 
(4.6) and taking the limit as R — > 0, one finds 



lim ( — 

R^o V E 



lim 

R^O 



(Co - po) + (Ci - Pl)R + (C2 ~ P2)R 2 + (Cs ~ P3)^ 3 + - 

(Co - po) + (Ci - \pi)R + (C 2 - Ip2)R 2 + (Cs - f ps)^ 3 + 



In order to avoid a divergence of (4.27) it is necessary to impose the condition 

Co = Po- 



(4.27) 



(4.28) 



This condition implies that no pointlike dipole can exist at R = 0, if one requires a finite 
E' /E value there. Whether a divergent E' /E at the origin is physically acceptable or not is 
yet to be seen. However, since we are interested in numerical calculations, we will require a 
finite origin value for E' / E. Applying (4.28) to (4.27), gives 



lim ( — 

R^O \ E 



lim 

R-^0 



(Ci - pi) + (C2 - P2)R + (Cs - ps)R 2 + (C 4 - p±)R 3 + ... 



(4.29) 



(Ci - {pi)R + (C 2 - § P 2)R 2 + (Cs - § ps)R 3 + ... 

In order to avoid a divergence of (4.29) it is necessary to impose the additional condition 

Ci = Pi, (4.30) 

which gives 



lim ( — 

R-+o \ E 



lim 

R^o 



(C 2 - P2) + (Cs - ps)R + (C4 - pi)R 2 + (Cs - ps)R 3 + - 

_ \P1 + (C2 - IP2)R + (Cs - lp 3 )R + (C4 - f P4)R 3 + ... 
C2 - P2 



Pi 



Pi 7^0. 



(4.31) 
(4.32) 



So with the conditions (4.28), (4.30) and p\ 7^ satisfied, E'/E will have a finite origin 
value given by (4.32). However, (4.30) also implies that such an arrangement would have 
a non-smooth central density profile, or 'cusp', since p\ ^ 0. In the case of smooth central 
density one has 

Pi = 0, (4.33) 
and so in order to avoid the divergence of (4.31), one must impose the further condition 

C2 = P2- (4.34) 
Applying (4.33) and (4.34) to (4.31) yields 



lim ( — 

>o\E 



lim 

r-»0 

5 
2 



(Cs - ps) + (C4 - Pi)R + (Cs - P 5 )R 2 + (Ce - Ps)R 3 + - 

§P2 + (Cs " IP3)R + (C4 - f PA)R 3 + (Cs - |P5)« 4 + ... 



C3 - P3 
P2 



P2^0. 



(4.35) 
(4.36) 



So, if (4.28), (4.30), (4.33) (4.34) and pi 7^ are all satisfied, the model has a smooth central 
density profile and the E'/E value at the origin is finite, and given by (4.36). Applying the 
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same reasoning as above to higher orders, one finds that the limiting value, in general, is 
given by 



lim '. 

r->0 \ E 



n + 3 ( Cn+l - Pn+l 



n \ p n 



Cn = Pn^0 (4.37) 



where n is the power in R of the first non-zero term of pit in (4.22). This implies that for a 
finite E'/E origin value, one requires p m in = Plt up to n th order in R. Also, for a non-zero 
origin limit, Cn+i 7^ Pn+i- In order to avoid having to choose p m in in such a way, we instead 
define it in terms of piT and a 'deviation function', as follows. 

PminiR) = Plt(R) [1 - p(R)] (4.38) 

where p satisfies 

< p(R) < 1 p{0) = (4.39) 

and 

^| fl =o = i = l..n (4.40) 

Choosing any function /i in this way is sufficient to ensure that the origin value of E'/E is 
finite. If, in addition to (4.39) and (4.40), one has 

lid* ' < 4 - 41 > 

the origin value of E'/E will be given by (4.37). Clearly, if (4.41) is not satisfied, E'/E will 
have an origin value of zero 10 . 

4.3 Expressing p max in terms of p min and p LT 

It will be useful to have an expression for the maximum density profile, p ma x, in terms of the 
density profiles, p m in and plt-, that does not involve E'/E. Rearranging (4.5) such that p is 
the subject, and separating the right hand side into two terms, one finds 

^{pLT-pAv) . . 

P = PLT + K _ & (4.42) 

R E 

Now, evaluating (4.42) at the density maximum and exploiting the anti-symmetrical property 
of E'/E (3.24), leads to 

~E \Pmi 



§\ Pmm (pLT ~ PAY) . . 

plt blT1E\ ( 3) 

R ' E \Pmin 

Substituting for E' /E\ Pmin from (4.5) with p = p m i n , and simplifying, one finds 

(PLT ~ Pmin){pLT ~ PAv) ... 

Pmax = PLT + 7 r— 7 T- (4.44) 

{PAV — Pmin) + {PLT ~ Pmin) 

So, by specifying the minimum and. LT density profiles, Pmin 

(r) and Plt{t\ °ne fixes the 

maximum density profile, Pmax(r)- 



10 It is worth noting that the origin value of E'/E is a coordinate dependant quantity (since it contains a 
prime), and the treatment above is only valid for the coordinate choice r = R. In our original formulation 
we made the coordinate choice r — M, which turned out to be far more restrictive on the allowable density 
profiles which produce a finite origin value of E'/E. 
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4.4 Evolution considerations 

4.4.1 LT density 

The LT-density, plt, will be specified on the initial and final time slices, as was done in [19]. 
When considering the time evolution of the model, the value of pet on intermediate time 
slices is calculated according to (A. 26) and the evolution of a at the origin found from the 
applicable one of (A.23) (A.24) and (A.25). 

4.4.2 Reconstructing p m in from (E' / E) max during evolution 

When reconstructing the model evolution it will be convenient to have an expression for Pmin 
in terms of (E' /E) max , instead of being in terms of {E' / E) Pmin . This is because (E' /E) max 
is constant in time. From the analysis of §4.1.2, we have 

5) = * (S) < 4 - 45) 

' Pmin \ ' max 

where 

X = -sign[R'(p LT - p AV )] (4.46) 
Hence, setting p = p m { n in (4.42), and substituting in (4.45), gives 

x(ir) (plt-pav) 

Pmin = PLT + ( EL\ ' ^ 

R X { E ), 



' max 



Similarly, for p ma x, one can write 

Mir) (PLT-PAV) 

_ „ \ J max (a aq\ 

Pmax - PLT 7^— T~E[\ ' ^ 

R + X V E I max 

Therefore, having calculated the profiles (E' /E) max (M), R(t, M), R'(t,M) and p LT (t,M) 
one may reconstruct the evolution of p m i n (t,M) and p max (t, M). 

4.4.3 Approximating R1/R2 

Since we choose to specify, at the final time, the density profiles and dipole parameters as 
a function of i?2, it will be useful to know the approximate range of R± which contains the 
same mass, prior to integrating the final density profile. This will allow one to specify the 
fluctuation length scale of the initial LT density profile, relative to the same set of worldlines. 
Thus, an approximate expression for R1/R2 is needed. Rearranging (3.31) and evaluating at 
the initial and final time, and dividing the two, leads to 



PAV 



Assuming the simulated density profiles are asymptotically uniform at large R, approaching 
some constant value, say pbg, an d if all the fluctuations are contained within the simulation, 
then we can make the rather crude approximation at the boundary of the simulation 

PAV ~ Pbg, (4.50) 
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and thus (4.49) gives 



e> / \ 1/3 

Ri„(PbqA . (4 . 51) 



R2 \Pbg,i 

Alternatively, recalling that for a flat FLRW evolution during the radiation dominated era 
the scale factor grows like a oc i 2 / 3 . If we assume that our simulation is not too different from 
parabolic evolution at the extremity of the simulation then we can approximate the growth 
rate of the areal radius to that of the scale factor from above. Hence, one can write 

(4.52) 



R2 \*2/ 

So, either (4.51) or (4.52) can be used to approximate the range of R±, given R2. 
4.5 The algorithm 

The above analysis has equipped us with the tools to construct the Szekeres model that 
evolves between two time slices, given initial and final density data. The procedure is broken 
into two parts. Firstly, obtaining the six arbitrary functions which define the metric from 
the specified initial and final data. And secondly, calculating the model evolution from those 
arbitrary functions. To begin, one must specify the following quantities: 

• Plt,i(Ri)'- The LT-density profile at t = t\ 

• Plt,2(R2)'- The LT-density profile at t = ti 

• ^{R-i): The 'deviation function' at t = ti 

• ® Pmi n {R2)'- The density minimum orientation angle, 6, at t = £2 

• 4>p min (R2) : The density minimum orientation angle, 0, at t = £2 

The approximate range of plt,i is related to that of plt, 2 by (4.51) or (4.52). Now, the 
metric functions M(i?2), /(-R2), *b(-R2), £(-^2)1 P{R2) and Q{R2) are obtained as follows: 

Evaluate the integral (2.15) (with p = plt) at the initial and final times, to find M(R\) and 
M(i?2)- Interpolation then gives R\{M) and i?2(M), for a set of worldlines. Now, for each 
particle worldline (at each M value): 

• Calculate a\{M) and a2(M) using (A.l). In the case where M m j n = 0, the origin value, 
a(0), is given by (A. 22) (with p = plt)- 

• Evaluate the inequalities (A. 2) (A. 6), (A. 10), (A. 14) (A. 18), and hence determine the 
evolution type of that worldline. 

• Determine the value of x which solves the relevant ^{x) = equation, (A. 5) (A. 9) 
(A. 13) (A. 17) (A. 21). A bisection method is good for this - refer to Table 1 for the 
initial rage over which to bisect. 

• Calculate f(M) using the one of (A.3) (A.7) (A. 11) (A.15) (A.19) which is applicable 
to the evolution type. 
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• Calculate t b {M) using the one of (A.4) (A.8) (A.12) (A.16) (A.20) which is applicable 
to the evolution type. 

• Calculate E'/E\ max (M) on the final time slice using (4.6) with (4.38) and (3.31). At 
the origin pAv{Q) = Plt(0), and E'/E\ max (0) is given by (4.37) if (4.41) is satisfied. If 
(4.41) is not satisfied, the origin value of E' / E\ max is zero. 

• Determine the sign of R'(plt-PAv) at final time, and hence calculate Sill yJrnax j S1H (pmax j 
cos9 max and cos4> max using either (4.7) or (4.8). 

• Calculate S, P and Q using (4.19), (4.20) and (4.21). 

Having determined the metric functions one can now check for regular origins, spatial maxima 
or minima and shell crossings (see §6.1 in [38]). Going on to reconstruct the model evolution, 
on a set of intermediate time slices, one calculates 

• a{M) using the one of (A. 23) (A. 24) (A. 25) which is applicable to the evolution type. 

• plt(M) using (A.26). 

• R(M) using (A.l). 

• R'(M) using (2.12). 

• pav(M) using (3.31) 

• x( M ) usin g (4.46). 

• pmin(M) and p max {M) using (4.47) and (4.48). 

5 Numerical simulations 
5.1 Implementation 

A script was written in MATLAB to implement the algorithm outlined in §4.5. Some practical 
aspects of the implementation are outlined below. 

• Evaluation of the integral in (2.15), to obtain M(R\) and M(i?2), is computed using 
the adaptive Gauss-Kronrod quadrature package, quadgk. The inverses, R\{M) and 
i?2(M), defined over approximately the same mass rage, are then found. The function 
interpl is used to interpolate the values of Ri at the same M values for which i?2 is 
defined. 

• The quantities a\ and pav,i are best calculated using the values of R and M given 
by the integration of (2.15), prior to interpolation. Once cq and pav,i are known, the 
function interpl is the used to interpolate values corresponding to the same M values 
for which R2 is defined. 

• The function m-file solve-phi_x.m was written to solve the various ip(x) = equations 
by the bisection method. 
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• Evaluation of the integral in (4.19) is computed using the ODE solver ode45, while 
the integrals in (4.20) and (4.21) are computed using the adaptive Simpson quadrature 
package quad. The function m-file myinerpfun.m was written so that the intergrands 
in (4.19) (4.20) and (4.21) are continuously defined, and thus amenable to integration 
with ode4-5 and quad. 

• For non-parabolic world lines, the evolution of R is not given explicitly in terms of t, 
but instead it is parametrized by rj. This necessitates the interpolation of R on constant 
time slices. Along constant-M slices an N x 1 vector of rj- values is created ranging from 
7]i to r/2 11 , allowing R{jf) and the corresponding t{rj) to be calculated. A spline cure 
is fitted to R{t) at each M value using spline, allowing interpolation of R at constant 
time slices. 

• Derivatives, such as R'(t,M), R(t,M) and a a/(*, M), are computed using the package 
gradient, which uses a 3-point finite difference method. 

5.2 Model profiles 

Since the aims of this paper are to develop a model construction procedure and implement it 
in software, investigating elaborate models using this software is thus outside our scope. Our 
choice of model profiles was therefore predominantly motivated by convenience, the intention 
being to sufficiently demonstrate the functioning of the software. We use geometric units, 
in which c = 1 = G and fix the scale freedom of GR by considering 1O 15 M0 to be one mass 
unit, as in [19]. The corresponding geometric length and time units are related by 

1M G = 10 15 Mq, 

=> ILq = Mq^k = 48 pc, 
cr 

^ir G = M G § = 156yr. (5.1) 

In all of the simulations we choose the final time to be, approximately, the current age of the 
Universe, and the initial time to be, approximately, the time of recombination. That is 

h = 100 k yr = 641 T G , (5.2) 
t 2 = lOGyr = 6.41 x 10 7 T G . (5.3) 

Similarly, background densities are given by 

PBG,l = 8 x 10~ 17 kg/m 3 = 1.3 x 10~ 7 M G /L G , 

Pbg,2 = 8 x 10~ 27 kg/m 3 = 1.3 x 10~ 17 M G /L 3 G . (5.4) 

In all cases, the range of R2 over which we simulate the chosen profiles is 10 6 Lq, which is 
equivalent to 48 Mpc. 

As developed above, the models will be specified in the following way. On the initial 
time slice, the "LT-density" is given. This shows the variation of the "equatorial" density (the 
density where E' / E = 0) with areal radius at that time. A similar profile is given for the final 
time. The strength of the density dipole is given at the final time via (4.38), as a deviation 

n For regions that are recollapsing at t = ti, the final phase is given by 772 = 2-7T — arccos(l — a^x) 
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from the "LT-density" , while two more profiles specify the (f> orientation of this minimum. 
Thus as a rough guide, the locus of shell-minimum density follows (R2, 6p min -, 0p m i„)- 12 

5.2.1 Run #1 

We first investigate the simple case where the dipole orientation is unchanging and E'/E\ max (0) 
0, and thus choose the 'deviation function' such that it does not satisfy (4.41). The chosen 
profiles are 

1.00003 + (8 x 10- 5 )fl 2 



PLT,l(Rl) 


= Pbg,i 




= PBG,2 




= 10~ 23 




7T 

" 4 




7T 

= 4' 



1 + (8 x 10~ 5 )i? 3 
10 5 + 10" 10 i? 2 + 10" 17 i? 3 
10 + 10" 10 ^ 2 + io- 17 i? 3 
J R 4 exp(-(4x 10- 24 )fi 4 ) 



(5.5) 
5.2.2 Run #2 

Next, we investigate the case of i?-dependant dipole angles and where E' /E\ max (0) / 0, and 
thus choose the 'deviation function' such that it does satisfy (4.41). The chosen profiles are 

.1.00003 + (10- 5 )i? 3 

PLT,l{Rl) = PBG,1 



PLT,2(R2) = PBG, 



1 + (10~ 5 ).R 3 
(5 x 10 4 ) + 10- 10 i? 2 + W- 17 R 3 



10 + 10~ 10 i? 2 + io- 17 i? 3 
p(R 2 ) = 10~ 17 i? 3 exp (-(4 x 10~ 18 )i? 3 ) 



<t> Pmin {R2) = (2vr x 10- 6 )i?, (5.6) 

5.2.3 Run #3 

Though the profiles are superficially similar to Run ^2, this case has a strongly varying late 
time Pmax profile. The chosen profiles are 

.1.00003 + (1 x lO- 5 )^ 3 

PLT,l\R\) = PBG,1 



PLT,2(Rz) = PBG 



1 + (1 x 10~ 5 )i? 3 
1 + (2 x 10~ W )R 2 + (2.9 x 10~ 17 )i? 3 



10 + (2 x 10~ 10 )i? 2 + (2.9 x 10~ 17 ).R 3 
p{R 2 ) = (2.2 x 10- 18 )i? 3 exp(-(1.2 x 10~ 17 )i? 3 ) 



o Pmin (R2) = J + (^xio- 6 )i? 



Ji? 2 ) = (27rxlO- 6 ) J R, (5.7) 



Pmi' 



12 A detailed visualisation would require a careful treatment of the non-concentricity of the shells, and would 
be complicated by the fact that the space is curved in general, while any plot would be in flat space. 
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5.3 Results and discussion 

5.3.1 Run #1 

The density profiles plt, Pmin, Pmax and pav at the initial and final times are shown in 
Figures 1 and 2, respectively. Both have an over-density at the origin, resulting in pav > Plt 
for all (t, r), and thus the relationship between E'/E\ max and E'/E\ Pmin is constant. At the 
final time, the large deviation in p m in away from plt (which was specified via the function 
p) causes little deviation in p ma x away from plt- The large deviation in p m in at the final 
time is sourced by a much smaller deviation from the initial time. The (8, 4>) orientation of 
the density dipole is constant. The function E' /E\ max shown in Figure 3 has an origin value 
of zero, as expected from our choice of p. Figure 4 shows the LT arbitrary functions, M, f 
and tb- Since there are no regions where plt = 0, we see that M in the top pane of Figure 
4 is always increasing. The LT energy function in the centre pane has / < 0, and thus all 
worldlines are elliptical. The bottom pane shows the bang time, which is decreasing close 
to the origin and then becomes increasing further away. This will produce shell crossings, 
however, given the gradual slope we expect these to occur very soon after the bang, well 
before t±, when we anticipate our model will not be valid anyway, due to its dust equation of 
state. The Szekeres arbitrary functions, S, P and Q, are shown in Figure 5. Surfaces showing 
the full time-evolution of plt, Pmin and p ma x are shown in Figures 6 7 and 8, respectively. 

5.3.2 Run #2 

The density profiles plt, Pmin, Pmax and pav at the initial and final times are shown in 
Figures 9 and 10, respectively. As in the case of Run $4, both have an over-density at the 
origin, resulting in pav > PLT for all (t,r), and thus the relationship between E'/E\ max and 
E'/E\ Pmin is constant. In this case the (9,(f>) orientation of the density dipole varies as R 
increases, describing one circuit of a helix over the range of the model. Again, the large 
deviation of p m in away from plt, at the final time, causes little deviation of p ma x away from 
Plt- The function E' /E\ max , shown by the broken line in Figure 3, approaches a finite origin 
value, as expected from our choice of p. The broken lines in Figure 4 shows the LT arbitrary 
functions, M, / and tb, which are similar to those from Run $ 1. The Szekeres arbitrary 
functions, S, P and Q, are shown by the broken lines in Figure 5. The effect of varying the 
dipole orientation angles is evident in all of them. Surfaces showing the full time evolution 
of plt, Pmin and p ma x are qualitatively very similar to those from Run #1 (see Figures 6, 7 
and 8), and are thus omitted here. 

5.3.3 Run #3 

The density profiles plt, Pmin, Pmax and pav at the initial and final times are shown in 
Figures 11 and 12, respectively. Unlike the previous two cases, the initial over-density evolves 
into an under-density at the final time, causing the quantity (plt~ Pav) to change sign during 
the course of the evolution. Thus, the relationship between E' /E\ max and E' /E\ Pmin is not 
constant for all (t, r) - a 'flip' takes place, occurring on a non-linear locus at or near t\, and 
covering the whole range of r within 1700 years. Though the density dipole orientation is 
as for Run #2, the density dipole orientation is reversed at the initial time. At the final 
time, the 'deviation function ' produces a large difference between plt and p m ax, but little 
between plt and p m in- This effect seems to be depend on the sign of (plt — Pav)- If 
R'(plt — Pav) > 0, then the deviation in maximum density is enhanced as compared to the 
minimum density, and vice versa. Furthermore, the location of the peak in p m ax corresponds 
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Density Profiles at t = t\ 




1 .00000 - 



ggggg 1 1 1 1 1 1 1 1 

0.002 0.004 0.006 0.008 0.01 0.012 

R\ [Mpc] 

Figure 1. Initial Density Profiles from Run #1 - The minimum, maximum, internal average 
and LT density profiles at t = t\ from Run Jfl. The LT-density was specified and the rest were 
calculated. (The profiles extend out to Ri = 0.27) 

to the location of the minimum in p m i n , which, by definition, corresponds to the peak in the 
deviation function, fi. The function E' /E\ max , shown in Figure 13, approached a finite origin 
value, as expected from our choice of p. Figure 14 shows the LT arbitrary functions, M, 
f and tf,. The Szekeres arbitrary functions, S, P and Q, are shown in Figure 15. Surfaces 
showing the full time evolution of ptT, Pmin and p ma x are shown in Figures 16 17 and 18, 
respectively. 

6 Conclusions 

Modern cosmology, although having enjoyed a great number of successes in the last century, 
is still faced with many open questions (e.g. the apparent dimming of distant supernovae). 
This has given rise to many highly speculative theories receiving much of the spotlight in 
recent years. While the importance of such research cannot be discounted, we should not 
loose sight of the full implications of the best current theory, GR. The consequences of the 
non-linearity of the EFE's have not yet been fully explored in cosmology, and investigating 
exact solutions is an indispensable tool for doing this. In this respect, the Szekeres family of 
inhomogeneous solutions offers a wide range of possibilities for modelling cosmic structure. 
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Density Profiles at t = t? 
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Figure 2. Final Density Profiles from Run #1 - The minimum, maximum, internal average 
and LT density profiles at t = t% from Run #1. The minimum and LT-dcnsity were specified. The 
rest were calculated. 
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Figure 3. The Function E'/E\ max from Run #1 and Run #2 - The dipole function E'/E\ r 
from Run #1 (solid line) and Run #2 (dashed line). 
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Figure 4. Arbitrary Functions A/, / and % from Run #1 and Run #2 - The arbitrary 
functions which are common to both Szekeres and LT modes; M, f and tt,, from Run #1 (solid line) 
and Run #2 (broken line) 



When attempting to construct a realistic model of the universe, or part thereof, it is of great 
utility to do so from physical quantities or data more directly accessible to observation than 
theoretical metric functions. In other words, there's a need to be able to create models that 
incorporate observational constraints from two different times, such as recombination and 
recent times. 

We considered quasispherical Szekeres models, outlining a model construction procedure 
using given density data at some initial and final time. Szekeres models are thought to 
be a sequence of non-concentric mass shells, each with density dipole. Consequently, the 
procedure requires one to specify 'radial' profiles of the equatorial density at the initial and 
final time, as well as dipole parameters (encoding the intensity and orientation), which we 
choose to specify at the final time. Since, in Szekeres models, the dynamics of the areal 
radius is identical to that in LT models, the evolution of each shell can be determined by 
the equatorial density profiles, pLT(ti,R) and Plt&i R), as is the case in LT models. We 
used a minor modification of the LT model construction procedure of Krasinski &; Hellaby 
to determine the arbitrary functions which are common to both - M, f and %. The dipole 
variations over each shell are encoded in the 'deviation function', [J,(R), that defines the 
minimum density p m in on each R shell, and the orientation angles, Pmin (R) and <ftp min (R), 
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Figure 5. Arbitrary Functions S, P and Q from Run #1 and Run #2 - The arbitrary 
functions which are unique to Szckcres models; 5, P and Q, from Run #1 (solid line) and Run #2 
(broken line) 

of that minimum density. With knowledge of plt(R) and fi(R) one can determine p m in(R) 
and Pmax(R)- For the determination of the metric functions unique to Szekeres models, S, P 
and Q, representing degrees of freedom pertaining to the dipole, we derive a new result - exact 
analytic expressions in terms of the dipole parameters, Pmin , 4> Pmi „ and E'/E\ max (which is 
directly related to p m in and plt)- So, by specifying the profiles pLr(ti, R), Plt^i, R)-, m(^)> 
Pmin (R) and <f>p min (R), and following the algorithm given in §4.5, one can determine the 
six arbitrary metric functions which completely define the Szekeres model. Special attention 
was paid to the origin limit of the dipole parameter E'/E, investigating the conditions which 
ensure it is finite and non-zero there. We corrected a claim made in previous literature that 
a maximum in E'/E corresponds to a density minimum (see Equation 69 in [38]) by showing 
that the derivative of the density with respect to E'/E is not always negative, but rather, it 
depends on the sign of R'(piT — PAv)- We found that a maximum in E'/E corresponds to 
a density minimum if R'(plt — Pav) < 0. 

Using MATLAB, code was written to implement the procedure for determining these 
metric functions, as well as to simulate the model evolution. Since investigating elaborate 
models is beyond the scope of this work, we considered only three simple cases. All models 
spanned the time from recombination until the present, with the choice of initial density 
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Figure 6. Evolution of the LT-density from Run #1 - The evolution of the LT-density profile, 
Plt- Log scale 
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Figure 7. Density Dipole Minimum from Run #1 - The evolution of the density dipolc 
minimum, p m in- Log scale 
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Figure 8. Density Dipole Maximum from Run #1 - The evolution of the density dipole 
maximum, p m ax- Log scale 




Rr [Mpc] 



Figure 9. Initial Density Profiles from Run #2 - The minimum, maximum, internal average and 
LT density profiles at t = t\ from Run #2. Only LT-density was specified. The rest were calculated. 
(The profiles extend out to R\ = 0.27) 
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Density Profiles at t = t 2 
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Figure 10. Final Density Profiles from Run #2 - The minimum, maximum, internal average 
and LT density profiles at t — t-2. from Run $=2. The minimum and LT-density was specified. The 
rest were calculated. 

profile consistent with CMB anisotropies. We then chose different profiles and orientation 
angles at the final time for each of the models. In all cases, when reconstructing the model 
evolution, the calculated metric functions reproduced the specified initial and final density 
profiles. In Runs #1 and #2 the initial and final profiles both have an over-density at the 
origin, causing the sign of (plt~ PAv) to be negative for all (t, r), and hence E' /E\ max occurs 
at density minimum for all (t,r). In both cases, the deviation of p m in away from p^x was 
much greater than that in p m ax- In contrast, Run ^3 had an initial over-density evolving 
into an under-density at the final time, which caused (plt — PAv) to change sign during the 
course of the model evolution. Also, the deviation in p m in away from plt, at the final time, is 
much less than the deviation in p m ax ■ We speculate that the profile which shows the biggest 
deviation away from pur, either Pmin or Pmaxi seems to depend on the sign of (plt — PAv)- 
In the future, plotting the density on a section of geodesic spatial 2-surfaces will give 
one a better feel for the resulting density distribution. Since the constant- (t, r) mass shells 
are interpreted as being arranged non-concentrically, the distance from the origin to a given 
shell is not trivially related to areal radius, but instead, has (9, (f>) variations. Hence, plotting 
Pmin, Pmax and plt as a function of r, whatever one's choice for r may be, does not give 
one a sense of the density profile in any direction. In addition, the local (0,4>) coordinates, 
which describe the dipole orientation on each shell, are not necessarily parallel. In order to 
understand the relationship between the coordinates on each shell some parallel transport 
operation may need to be performed. It is also foreseeable that, in the future, the model 
construction algorithm could be extended to include velocity profiles, as was done in [20] for 
the LT model. Extending these results to the e / 1 Szekeres cases will require care as the 
interpretation of R and M is different from the quasi-spherical case [37] . 
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Density Profiles at t = t\ 




Figure 11. Initial Density Profiles from Run #3 - The minimum, maximum, internal average 
and LT density profiles at t = t\ from Run #3. Only LT-density was specified. The rest were 
calculated. 

A Details of the LT model construction 

Here we provide some important equations from the LT model construction procedure of [19], 
in which the evolution is determined from initial and final density profiles, though adapted 
to the coordinate choice (2.10). In particular, we list the cases that arise, their ranges of 
validity, and the equations to be solved in each case. This is a core part of the proposed 
Szekeres model construction procedure presented in §4. 

A.l Finding / and if, 

This material follows on from (2.15) of section 2.2. 
With the variables 13 

— Ri — IZJ (\ -n 

the nature of the LT model that evolves between the initial and final time slice at a given M 
is 

Hyperbolic (/ > 0) 

if 

t2-h<^(a 3 2 /2 -af), (A.2) 
then, the energy function is given by 

/ = xM 2/3 , (A.3) 



13 a and x have the advantage of being non-zero at the origin 
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Figure 12. Final Density Profiles from Run #3 - The minimum, maximum, internal average 
and LT density profiles at t = ti from Run #3. The minimum and the LT-density was specified at t\ 
and ti- The rest were calculated. 
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Figure 13. The Function E' / E\ max from Run #3 - The dipole function E'/E\ max from Run #3 
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Figure 14. Arbitrary Functions M, f and i& from Run #3 - The arbitrary functions which 
are common to both Szekeres and LT modes; M, / and tb 



and the bang time by 

— % 

where x solves 



y/ (1 + xai) 2 — 1 — arcosh(l + ajx) 



(A.4) 



= ipHx{x) =y (1 + a,2x) 2 — 1 — arcosh(l + a 2 x) 

— \/ (1 + aix) 2 — 1 + arcosh(l + a\x) — (t 2 — t\). 



,3/2 



(A.5) 



Near Parabolic (/ 0) 

if (t 2 —tx) is close to 



t 2 -h = y^(4 /2 -al /2 ), (A.6) 



3 

then, the energy function is given by 

/ = xM 2/3 , (A.7) 

and the bang time by 

3/2 ( ' 3 9 _ 



t 6 = ti - < /z ( 1 - -a^ + — a^V ) , (A.8) 
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Figure 15. Arbitrary Functions S, P and Q from Run #3 - The arbitrary functions which are 
unique to Szckcres models; S, P and Q 



where x solves 



= iP P {x) 



V2X 3 ? 2 



3/2 



9 



2„,2 



aJ I 1 a 2 x H tO,2 x 



20 224 



3/2 



1 ~ 20 aiX + 224°^ J — (*2 — *i) 



(A.9) 



Elliptic and still expanding at t 2 (/ < 0) 

if 



^2 2 . 2 

(a 2 /2) 3 / 2 7r - arccos(l - 2ai/a 2 ) + 2^/ai/a 2 - {ai/a 2 ) 2 > t 2 - h > — (a^ 2 - a^ 2 ) 



then, the energy function is given by 



and the bang time by 

tb — ^2 ^ 



/ = - X M 2/3 , 



arccos 



(1 - dix) - v 7 ! - (1 - a%x) 2 



(A.10) 
(A.ll) 

(A.12) 
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Figure 16. Evolution of the LT-density from Run #3 - The evolution of the LT-density profile, 
Plt- Log scale 
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Figure 17. Density Dipole Minimum from Run #3 - The evolution of the density dipolc 
minimum, p m in- Log scale 



- 35 - 



Maximum Density 



5 -30 




10 50 



t [yrs] 



R 2 [Mpc] 



Figure 18. Density Dipole Maximum from Run #3 - The evolution of the density dipolc 
maximum, p m ax- Log scale 



where x solves 



= ipx{x) = arccos(l — d2x) — yl — (T— a 2 x)' 



arccos 



(1 - ai x) + V 1 - (1-aix) 2 - (t 2 - h). 



,3/2 



(A.13) 



Elliptic and near maximum expansion at t 2 (/ < 0) 

if (t2 — ti) is close to 



ti — t\ = (a 2 /2) 3 / 2 7r - arccos(l - 2ai/a 2 ) + 1\J a\ja 2 - (ai/a 2 ) 2 
then, then energy function is given by 

/ = -xM 2 l\ 

and the bang time by 



tb = t\ — x 



t 2 - x 



-3/2 



-3/2 



where x solves 



arccos(l — a\x) — yl — (1 — a\x) 2 



2 3/2 

7T _ 2 3 / 2 (2 - asrr) 1 / 2 + —(2 - a 2 x) 



2 3/2 



3/2 



(A.14) 
(A.15) 



(A.16) 



= Vm(x) =vr - 2 3 / 2 (2 - asx) 1 / 2 + —(2 - a 2 x) 3 / 2 

— arccos(l — a\x) + y 7 aix(2 — a%x) — (t 2 — ii)a; 3 / 2 . (A.17) 
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Elliptic and recollapsing at t 2 (/ < 0) 

if 

h — h > ifl^l'ifl 2 7r - arccos(l - 2ai/a 2 ) + 2y l a\ja 2 - (ai/a 2 ) 2 



then, the energy function is given by 



and the bang time by 

t b = tt- x~ 3/2 
= t 2 - x- 3 / 2 

where x solves 



arccos(l — a±x) — yl — (T— a\x) 2 
it + arccos(— 1 + a 2 x) + \/l — (1 — 



= ipc(x) =7r + arccos(— 1 + a2.x) + \J\ — (1 — a^x) 2 

— arccos(l — aix) + \/l — (1 — aix) 2 — (i 2 — t±)x 

A. 2 Solving -0(a;) = 



3/2 



(A.18) 
(A.19) 



(A.20) 



(A21) 



The value of x which solves tp(x) = 0, at each M value, can be found numerically using the 
bisection method. The range in x over which to bisect and a good starting value for the first 
for each evolution type, are give in Table 1. 



Evolution Type 


Bisection Range 


HX 


pi I an— ai \2 
to~U > 


all other types 





Table 1. Bisection Info - Showing, for various evolution types, the range for the bisection method 
which is used to solve the various ip(x) = equations. 

Increasing numerical error in the exact expressions for ip(x) necessitates the use of series 
expansions, as the borderline cases are approached. Thus, 'fat' borderline equations are used 
when solving for x in regions near the PX or EM cases. The range over which these 'fat' 
borderline expressions are valid is assumed to be, approximately, the region where the ratio 
of the magnitude of the third order component of the series expansion, to that of the first, is 
less than 10~ 3 . 

A. 3 Limiting values at M = 

While several quantities have the value at the origin, the variables and x have finite 
limits as M — > 0. The origin value of a is [19] 



Oi(0) 



G 



«Pi(0) 



1/3 



(A.22) 



and the value of x, calculated by solving the relevant ip{x) = equation, comes out non-zero 
automatically when non-zero values of aj(0) are used. 
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A. 4 Reconstructing model evolution 

When reconstructing the evolution of the model, it is convenient to re-write the LT solutions 
in terms of a and x. The time evolution of the quantity a, for the various evolution types, is 
given below. 



Elliptic: 



Parabolic or close to it: 



1 ~ c o s V . . (t; - sin 77) . . 

a = " ' t ~ t b= x s/2 ( A - 23 ) 



[6(t - t b )f 



504000 1 

(A.24) 



Hyperbolic: 



a= co^-l t (smhg - V ) 

x x A l z 



In all cases the density is given by 



P 4ira 2 (a/3 + Ma, M ) ^ A ' 2 ^ 
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